Observing the dynamics of super-massive black hole binaries with Pulsar Timing Arrays 
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Pulsar Timing Arrays are a prime tool to study unexplored astrophysical regimes with gravitational 
waves. Here we show that the detection of gravitational radiation from individually resolvable super- 
massive black hole binary systems can yield direct information about the masses and spins of the black 
holes, provided that the gravitational-wave induced timing fluctuations both at the pulsar and at the Earth 
are detected. This in turn provides a map of the non-linear dynamics of the gravitational field and a new 
avenue to tackle open problems in astrophysics connected to the formation and evolution of super-massive 
black holes. We discuss the potential, the challenges and the limitations of these observations. 
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Introduction — Gravitational waves (GWs) provide a 
new means for studying black holes and addressing open 
questions in astrophysics and fundamental physics: from 
their formation, evolution and demographics, to the as- 
sembly history of galactic structures and the dynamical 
behaviour of gravitational fields in the strong non-linear 
regime. Specifically, GW observations through a network 
of radio pulsars used as ultra-stable clocks - Pulsar Tim- 
ing Arrays (PTAs) - represent the only direct obser- 
vational avenue for the study of supermassive black hole 
binary (SMBHB) systems in the ~ 10 8 - 1O 9 M mass 
range, with orbital periods between ~ 1 month and a few 
years, see ej>. raja and references therein. Ongoing obser- 
vations ]4j, |5j |8l |9|1 an d future instruments, e.g. the Square 
Kilometre Array jflfJll. ar e expected to yield the necessary 
timing precision fill H2I1 to observe the diffuse GW back- 
ground. This background is likely dominated by the inco- 
herent superposition of radiation from the cosmic popula- 
tion of massive black holes lfl3l - l2l1l and within it, we ex- 
pect a handful of sources that are sufficiently close, massive 
and high-frequency to be individually resolvable 1 22 - 30Tl. 



Massive black hole formation and evolution scenar- 
ios [31-13] predict the existence of a large number of 
SMBHBs. Furthermore, SMBHs are expected to be (pos- 
sibly rapidly) spinning [35, 3(3] . In fact the dynamics of 
such systems - which according to general relativity are 
entirely determined by the masses and spins of the black 
holes [37] - leave a direct imprint on the emitted gravita- 
tional waveforms. From these, one could measure SMBHB 
masses and their distribution, yielding new insights into the 
assembly of galaxies and the dynamical processes in galac- 
tic nuclei Il26ll . Moreover, measuring the magnitude and/or 
orientation of spins in SMBHBs would provide new infor- 
mation on the role of accretion processes ll38l - l42ll . Finally, 
detections of SMBHBs could allow us to probe general rel- 
ativistic effects in the non-linear regime in an astrophysical 
context not directly accessible by other means, see ll43ll 
and references therein. 

The observation of GWs with PTAs relies on the detec- 



tion of the small deviation induced by gravitational radi- 
ation in the times of arrival (TOAs) of radio pulses from 
millisecond pulsars that function as ultra-stable reference 
clocks. This deviation, called the residual, is the differ- 
ence between the expected (without GW contribution) and 
actual TOAs once all the other physical effects are taken 
into account. The imprint of GWs on the timing resid- 
uals is the result of how the propagation of radio waves 
is affected by the GW-induced space-time perturbations 
along the travel path. It is a linear combination of the 
GW perturbation at the time when the radiation transits 
at a pulsar, the so-called "pulsar term", and thenjvhen it 
passes at the radio receiver, the "Earth term" fl-dl. The 
two terms reflect the state of a GW source at two different 
times of its evolution separated by r = (1 + Cl - p) L p ~ 

3.3 x 10 3 (1 + ft • p) (L p /1 kpc) yr where fl and p are 
the unit vectors that identify the GW propagation direction 
and the pulsar sky location at a distance L p from the Earth, 
respectively, see e.g. ll24Tl . [We use geometrical units in 
which G = c = 1.] In a network (array) of pulsars all 
the perturbations at the Earth add coherently and therefore 
boost the signal-to-noise ratio (S/N) of the signal. Each 
pulsar term is at a slightly different frequency since the or- 
bital period of the binary evolves over the time r. 

Measuring the key physics of SMBHBs is ham- 
pered by the short (typically T = 10 yr) observation 
time compared to the typical orbital evolution timescale 
/// = 1.6 x 10 3 (7W/10 9 M Q )- 5 / 3 (//50nHz)- 8 / 3 yr 
of binaries that are still in the weak field adia- 
batic inspiral regime, with an orbital velocity v = 
0.12(M/10 9 M Q ) 1 / 3 (//50nHz) 1 / 3 JH. Here M = 
mi + m2, n = mxm<ilM and M. = M 2 / 5 /i 3 / 5 are the to- 
tal, reduced and chirp mass, respectively, of a binary with 
component masses mi j2 , and / is the GW emission fre- 
quency at the leading quadrupole order. The chirp mass 
determines the frequency evolution at the leading Newto- 
nian order. In the post-Newtonian (pN) expansion of the 
binary evolution 114511 in terms of v <C 1, the second mass 
parameter enters at p 1 N order (0(v 2 ) correction); spins 
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contribute at p 15 N order and above (0(v 3 )) causing the or- 
bital plane to precess through spin-orbit coupling, at lead- 
ing order. These contributions are therefore seemingly out 
of observational reach. 

The GW effect at the pulsar - the pulsar term - may 
be detectable in future surveys, and for selected pulsars 
their distance could be determined to sub-parsec preci- 



in the TOA residuals. At the leading Newtonian order, A4 
drives the frequency and therefore the phase <$(i) evolu- 
tion, with the second independent mass parameter entering 
from the p : N onwards. SMBHs are believed to be rapidly 
spinning, and the spins are responsible for three distinctive 
imprints in the waveform: (i) they alter the phase evolu- 
tion through spin-orbit coupling and spin-spin coupling at 



sion 1271 1461 14711 . If this is indeed the case, it opens the 
opportunity to coherently connect the signal observed at 
the Earth and at pulsars, therefore providing snapshots of 
the binary evolution over ~ 10 3 yr. These observations 
would drastically change the ability to infer SMBHB dy- 
namics, and study the relevant astrophysical process and 
fundamental physics. 

In this Letter we show that for SMBHBs at the high end 
of the mass and frequency spectrum observable by PTAs, 
say rrii 2 — 10 9 M Q and / = 10 -7 Hz, the observations 
of a source still in the weak-field regime become sensi- 
tive to post-Newtonian contributions up to p 15 N, includ- 
ing spin-orbit effects, if both the pulsar and Earth term can 
be detected. This in principle enables the measurement of 
the two mass parameters and a combination of the spin's 
magnitude and relative orientation. We also show that the 
Earth-term only can still be sensitive to spin-orbit coupling 
due to geometrical effects produced by precession. We dis- 
cuss the key factors that enable these measurements, and 
future observational prospects and limitations. 

Signals from SMBHBs — Consider a radio pulsar emit- 
ting radio pulses at frequency in the source rest-frame. 
GWs modify the rate at which the radio signals are re- 
ceived at the Earth UHll], inducing a relative frequency 
shift 8v(t)/v a = hit - r) - h(t), where h{t) is the GW 
strain. The quantities that are actually produced at the end 
of the data reduction process of a PTA are the timing resid- 
uals, J dt' 5u(t')/u , although without loss of generality, 
we will base the discussion on hit). The perturbation in- 
duced by GWs is repeated twice, and carries information 
about the source at time t, the "Earth term", and at past 
time t — t, the "pulsar term". 

We model the radiation from a SMBHB using the so- 
called restricted pN approximation, in which pN correc- 
tions are included only in the phase and the amplitude is 
retained at the leading Newtonian order, but we include 
the leading order modulation effects produced by spin-orbit 
coupling. The strain is given by 

h(t) = -A gw (t)A p {t) cos[$(t) + <p p (t) + tp T (t)] , (1) 

where A gw (t) = 2[7rf{t)] 2/3 M 5/3 /D is the Newtonian 
order GW amplitude, <J?(i) is the GW phase, see e.g. 
Eq. (232, 234) in JH and Eq. (8.4) in £1, and D is the 
distance to the GW source. A p (t) and ip p (t) are the time- 
dependent polarisation amplitude and phase and ip^ (t) is 
an additional phase term, analogous to Thomas precession, 
see Eq. (29) in JH]. 

The physical parameters leave different observational 
signatures in the GW strain h{t) and are therefore found 



p N and p N order, respectively 15011 . (ii) they cause the 
orbital plane to precess due to (at lowest order) spin-orbit 
coupling and therefore induce amplitude and phase modu- 
lations in the waveform through A p {t) and <£ p {t)\ and (hi) 
through orbital precession they introduce an additional sec- 
ular contribution </?t(£) to the waveform phase. Astrophys- 
ically we expect PTAs to detect SMBHBs of comparable 
component masses l24Tl . We therefore model the spin-orbit 
precession using the simple precession approximation ll49ll . 
which formally applies when rrii = m 2 , or when one of 
the two spins is negligible with respect to the other. Let 
Si. 2 and L be the black holes' spins and the orbital angu- 
lar momentum, respectively. Then both S = Si + S 2 and 
L precess around the (essentially) constant direction of the 
total angular momentum, J = S + L, at the same rate 
da/dt = it 2 (2 + 3m 2 /(2mi)) (|L + S|)/ 2 (t)/M S, 
where a is the precession angle, while preservin g th e angle 
of the precession cone, A^, see Fig. 4 of Ref. 14911 . This 
approximation is adequate to conceptually explore these ef- 
fects, however in the case of real observations, one will 
need to consider the exact expressions ||5lh . 

The detection and particularly the measurement of the 
aforementioned parameters relies on coherently matching 
the signal with a template that faithfully reproduces its am- 
plitude and, importantly, its phase evolution. We therefore 
consider the contribution to the total number of wave cycles 
a proxy for the significance of a specific parameter. Indi- 
vidual terms that contribute ~ 1 GW cycle or more mean 
that the effect is in principle detectable, hence one can infer 
information about the associated parameter(s). We show 
that information about the parameters can only be inferred 
for SMBHBs at the high end of the mass spectrum and PTA 
observational frequency range. Having a sufficiently high- 
mass and high-frequency GW source is also essential to en- 
sure sufficient frequency evolution over the time r, so that 
the Earth and pulsar term are clearly separated in frequency 
space cf. Table U We therefore take fiducial source param- 
eters of mi = m 2 = 10 9 Mq, frequency at the Earth at the 
beginning of the observation f E = 10~ 7 Hz and an obser- 
vational time T = 10 years to illustrate the main results. 
We provide scaling relations as a function of the relevant 
quantities, allowing the reader to rescale the results for dif- 
ferent astrophysical and/or observational values. 

Observations using the Earth-term only — We start by 
considering analyses that rely only on the Earth-term con- 
tribution to the residuals, as done in Ref. (25, 52]. The 
case of a coherent analysis based both on the Earth- and 
pulsar-term, introduced in Ref. ]22j], is discussed later in 
this Letter. Table U shows that, in general, the frequency 
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change over 10 yr is small compared to the frequency bin 
width, 3.2(10 yr/T) nHz JH HI. The observed sig- 
nal is effectively monochromatic, making the dynamics 
of the system impossible to infer. However, the presence 
of spins affects the waveform not only through the GW 
phase evolution, but also via the modulations of A p (t) 
and <p p (t) that are periodic over the precession period, 
and also introduces the secular contribution (pr(t). For 
m 12 = 10 9 M Q and f E = 10~ 7 Hz the orbital angu- 
lar momentum precesses by Aa = 2 rad (for dimension- 
less spin parameter a = S/M 2 = 0.1) and Aa = 3 rad 
(for a = 0.98), and therefore the additional modulation ef- 
fect on A p (t) and <p p (t) is small, and likely undetectable. 
However, the overall change of ^(t) over 10 yrs could 
be appreciable: the average contribution for each preces- 
sion cycle of this additional phase term is (Aip?) = 4ir 
or 4-7T (1 — cos depending on whether Cl lies inside 
or outside the precession cone, respectively l49h . If Q 
lies inside the precession cone, and given that the observa- 
tion will cover between a third and a half of a full preces- 
sion cycle, then (Atp-r) ~ tt, which could surely indicate 
the presence of spins. On the other hand, the precession 
cone will be small in general since \S/L\ ~ av{M/jj) ~ 
0.1 a (M//i) (M/10 9 Mq) 1 / 3 (//lOOnHz) 1 / 3 , therefore 
the likelihood of ft lying inside the precession cone 
is small, assuming an isotropic distribution and ori- 
entation of sources. In this case the Thomas pre- 
cession contribution (per precession cycle) is sup- 
pressed by a factor (1 — cosA^) ~ ^1/2 ~ 
5 x lO- 3 a 2 (M/fj,) 2 (M/W 9 M ) 2 / 3 (//lOOnHz) 2 / 3 , 
which will produce a negligible contribution A(p T (t) <C 1. 
However unlikely, spins may still introduce observable ef- 
fects that need to be taken into account. 

Measuring SMBHB evolution using the Earth and pulsar 
term — With more sensitive observations and the increasing 
possibility of precisely determining L p see e.g. ll47ll . the 
prospect of also observing the contribution from the pulsar- 
term from one or more pulsars becomes more realistic. We 
show below that if at least one of the pulsar terms can be 
observed together with the Earth-term, this opens opportu- 
nities to study the dynamical evolution of SMBHBs and, 
in principle, to measure their masses and spins. This is a 
straightforward consequence of the fact that PTAs become 
sensitive to ~ 10 3 yr of SMBHB evolutionary history, in 
"snippets" of length T L p that can be coherently con- 
catenated. 

The signal from each pulsar term will be at a S/N which 
is significantly smaller than the Earth-term by a factor ~ 
y/Np, where N p is the number of pulsars that effectively 
contribute to the S/N of the array. For example, if the Earth- 
term yields an S/N of ~ 36a/ N p /20, then each individual 
pulsar term would give an S/N ~ 8. The possibility of co- 
herently connecting the Earth-term signal with each pulsar 
term becomes therefore a question of S/N, prior informa- 
tion about the pulsar-Earth baseline and how accurately the 



SMBHB location in the sky can be reconstructed, as part 
of a "global fit", e.g. lf27ll . Assuming for simplicity that the 
uncertainties on L p and Cl are uncorrected, this requires 
that the distance to the pulsar and the location of the GW 
source are known with errors ^ 0.01(100 nHz //) pc and 
£ 3(100 nHz//)(l kpc/Lp) arcsec, resp ectively. These 
are very stringent constraints ll24l l28i l47ll . and a detailed 
analysis is needed in order to assess the feasibility of reach- 
ing this precision. Clearly if an electromagnetic counter- 
part to the GW source were to be found i53l l54ll . it would 
enable the identification of the source location in the sky, 
making the latter constraint unnecessary. 

We can now consider the contribution from the differ- 
ent terms in the pN expansion to the total number of cy- 
cles in observations that cover the GW source evolution 
over the time r that are encoded in the simultaneous anal- 
ysis of the Earth and pulsar terms. The results are sum- 
marised in Table U for selected values of m^, and f E 
and for a fiducial value r = 1 kpc. The wavecycle con- 
tributions from the spin-orbit parameter are normalised to 

P = (V 12 ) ELi [H3(m l /M) 2 + 75??] L • S u which 
has a maximum value of 7.8. Contributions from the p 2 N 
order spin-spin terms are negligible. The results clearly 
show that despite the fact that the source is in the weak 
field regime the extended Earth-pulsar baseline requires the 
p 15 N, and in some rare cases the p 2 N contribution, to accu- 
rately (i.e. within ~ 1 GW cycle) reproduce the full phase 
evolution. 

For mi 2 = 10 9 M Q and f E = 10~ 7 Hz there is a total 
of 4305 GW cycles over a 1 kpc light travel time evolu- 
tion, with the majority (4267) accounted for by the lead- 
ing order Newtonian term, providing information about the 
chirp mass, and tens of cycles due to the p J N and p 15 N 
terms (77 and 45, respectively), that provide information 
about a second independent mass parameter. Spins con- 
tribute to phasing at p 1 5 N with ~ 3/3 cycles. Therefore 
their total contribution is smaller than the p 15 N mass con- 
tribution by a factor between a few and ~ 10 . The addi- 
tional Thomas precession phase contribution may become 
comparable to the p J N mass contribution in some cases. In 
fact, for a = 0.1(0.98) the binary undergoes 24 (34) pre- 
cession cycles. This corresponds to a total Thomas preces- 
sion phase contribution of 306 (426) rad if f2 lies outside 
the precession cone. 

The modulations of A p (t) and <p p (t) are characterised 
by a small Al, because for most of the inspiral S <C L, 
and are likely to leave a smaller imprint on the waveform 
than those discussed so far. We can indeed estimate the 
importance of this effect for the most favourable parameter 
combinations. The value of (p p (t) oscillates over time with 
an amplitude which depends on the time to coalescence, S, 
L, Ct and p. We choose the orientation of S such that X L is 
maximised, and we vary Ct and p, each of which is drawn 
from a uniform distribution on the two-sphere. 

In Figure 1 we show that for rapidly spinning (a = 0.98) 
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mi(M ) m 2 (M Q ) / B (nHz) 


(v/c x 10" 


~ 2 ) timespan A/ (nHz) 


Total 


Newtonian 


p x N 


p 15 N spin-orbit//3 


p 2 N 


100 

10 9 10 9 


14.6 


10 yr 


3.22 


32.1 


31.7 


0.9 


-0.7 


0.06 


0.04 


9.6 


-1 kpc 


71.2 


4305.1 


4267.8 


77.3 


-45.8 


3.6 


2.2 


50 


11.6 


10 yr 


0.24 


15.8 


15.7 


0.3 


-0.2 


0.01 


< 0.01 


9.4 


-1 kpc 


23.1 


3533.1 


3504.8 


53.5 


-28.7 


2.3 


1.2 


100 

10 8 10 s 


6.8 


10 yr 


0.07 


31.6 


31.4 


0.2 


-0.07 


< 0.01 


< 0.01 


6.4 


-1 kpc 


15.8 


9396.3 


9355.7 


58.3 


-19.9 


1.6 


0.5 


50 


5.4 


10 yr 


0.005 


15.8 


15.7 


0.06 


-0.02 


< 0.01 


< 0.01 


5.3 


-1 kpc 


1.62 


5061.4 


5045.8 


20.8 


-5.8 


0.5 


0.1 



TABLE I: Frequency change A/, total number of GW cycles and individual contributions from the leading order terms in the pN 
expansion over the two relevant time scales - a 10 yr period starting at the Earth and the time period L p (l + ft ■ p) between the Earth 
and pulsar term (hence the negative sign) - for selected values of mi 2 and f^. 



SMBHBs this effect could introduce modulations larger 
than 7r/2 in ip p (t) over 30% of the parameter space of 
possible ft and p geometries. The amplitude would cor- 
respondingly change over the same portion of the parame- 
ter space by at most 60% with respect to its unmodulated 
value. Since these effects are modulated, they will not be 
easily identifiable. 

Conclusions — We have established that the coherent ob- 
servation of both the Earth and pulsar term provides in- 
formation about the dynamical evolution of a GW source. 
The question now is whether they can be unambiguously 
identified. A rigorous analysis would require extensive 
simulations based on the actual analysis of synthetic data 
sets. We can however gain the key information with a 
much simpler order of magnitude calculation. The phase 
(or number of cycles) error scales as ~ 1 /(S /N). Assum- 
ing S/N ~ 40 means that the total number of wave cycles 
over the Earth-pulsar baseline can be determined with an 
error ~ 4300/40 ~ 100 wave cycles. This is compara- 
ble to the p 1 N contribution to the GW phase and, in very 
favourable circumstances, to the Thomas precession phase 
contribution, and larger by a factor of a few or more than 
all the other contributions. It may therefore be possible to 
measure the chirp mass and, say, the symmetric mass ratio 
of a SMBHB, and possibly a combination of the spin pa- 
rameters. Effects due to the p 15 N and higher phase terms 
are likely to remain unobservable, as well as amplitude 
and phase modulations. Correlations between the param- 
eters, in particular masses and spins, will further degrade 
the measurements. The details will depend on the actual 
S/N of the observations, the GW source parameters, and 
the accuracy with which the source location and the pul- 
sar distance can be determined. We plan to explore these 
issues in detail in a future study. 
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FIG. 1: The fraction of parameter space in ft and p for which 
the maximum excursion of <^ p over the time L p (l + ft ■ p) for 
L p = 1 kpc exceeds a certain value, shown on the horizontal axis. 
Several values of m 12 , a and /e are considered (see legend). 
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